<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Example 6.6: Comparison of worst-case robust, Tikhonov, and nominal least squares</title>
<link rel="canonical" href="/Users/mcgrant/Repos/CVX/examples/cvxbook/Ch06_approx_fitting/html/wcrobls.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Example 6.6: Comparison of worst-case robust, Tikhonov, and nominal least squares</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
Text output
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.4.2, Figure 6.16</span>
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX Argyris Zymnis - 11/27/05</span>
<span class="comment">% (a figure is generated)</span>
<span class="comment">%</span>
<span class="comment">% Consider the least-squares problem:</span>
<span class="comment">%       minimize ||(A0 + u1*A1 + u2*A2)x - b||_2</span>
<span class="comment">% where u = [u1 u2]' is an uncertain parameter and ||u||_2 &lt;= 1</span>
<span class="comment">% Three approximate solutions are found:</span>
<span class="comment">%   1- nominal optimal (i.e. letting u=0)</span>
<span class="comment">%   2- Tikhonov Regularized Solution:</span>
<span class="comment">%           minimize ||A0*x - b||_2 + delta*||x||_2</span>
<span class="comment">%      for some delta (in this case we set delta = 0.1)</span>
<span class="comment">%   3- worst-case robust approximation:</span>
<span class="comment">%           minimize sup{||u||_2 &lt;= 1} ||(A0 + u1*A1 + u2*A2)x - b||_2)</span>
<span class="comment">%      (reduces to solving an SDP, see pages 323-324 in the book)</span>

m = 50;
n = 20;
randn(<span class="string">'state'</span>,0);
rand(<span class="string">'state'</span>,0);

A0 = randn(m,n);
[U,S,V] = svd(A0);
S= diag(fliplr(logspace(-0.7,1,n)));
A0 = U(:,1:n)*S*V';
A1 = randn(m,n);  A1 = A1/norm(A1);
A2 = randn(m,n);  A2 = A2/norm(A2);

Aperb0 = [A1;A2];
p = 2;

b = U(:,1:n)*randn(n,1) + .1*randn(m,1);

<span class="comment">% we consider LS problems || (A0 + u1*A1 + u2*A2) x - b||</span>
<span class="comment">% where  ||u|| leq rho</span>

<span class="comment">% Nominal Solution</span>
xnom = A0\b;

<span class="comment">% Tikhonov Regularized Solution</span>
delta = .1;
xtych =  [A0; sqrt(delta)*eye(n)] \ [b; zeros(n,1)];

<span class="comment">% Robust Least Squares solution</span>
cvx_begin <span class="string">sdp</span> <span class="string">quiet</span>
    variables <span class="string">t</span> <span class="string">lambda</span> <span class="string">xrob(n)</span>
    minimize(t+lambda)
    subject <span class="string">to</span>
        [eye(m) A1*xrob A2*xrob A0*xrob-b; <span class="keyword">...</span>
         [A1*xrob A2*xrob]' lambda*eye(2) zeros(2,1); <span class="keyword">...</span>
         [A0*xrob-b]' zeros(1,2) t] &gt;= 0;
cvx_end

<span class="comment">% Generate Random Trials</span>
notrials=100000;
r = sqrt(rand(notrials,1));     <span class="comment">% random on [0,1] with pdf g(r) = 2r;</span>
theta = 2*pi*rand(notrials,1);  <span class="comment">% uniform on [0,2pi]</span>
v = [r.*cos(theta)  r.*sin(theta)];
ls_res = zeros(1,notrials);
rob2_res = zeros(1,notrials);
rob_res = zeros(1,notrials);
tych_res = zeros(1,notrials);

<span class="keyword">for</span> i =1:notrials

  A = A0 + v(i,1)*A1 + v(i,2)*A2;
  ls_res(i) = norm(A*xnom-b);
  rob_res(i) = norm(A*xrob-b);
  tych_res(i) = norm(A*xtych-b);

<span class="keyword">end</span>;


<span class="comment">% Plot histograms</span>
figure
<span class="comment">%subplot(211)</span>
[N1, hist1] = hist(ls_res,[min(ls_res):.1:max(ls_res)]);
freq1 = N1/notrials;
[N2, hist2] = hist(rob_res,hist1);
freq2 = N2/notrials;
[N3, hist3] = hist(tych_res,hist1);
freq3 = N3/notrials;



h = bar(hist3,freq3);
text(3, 0.07, <span class="string">'Tikhonov'</span>);
set(h,<span class="string">'FaceColor'</span>,0.90*[1 1 1]);
hold <span class="string">on</span>

h = bar(hist2,freq2);
text(4.2, 0.05, <span class="string">'Nominal'</span>);
set(h,<span class="string">'FaceColor'</span>,0.80*[1 1 1]);

h = bar(hist2,freq2);
set(h,<span class="string">'FaceColor'</span>,<span class="string">'none'</span>);
text(2.6, 0.2, <span class="string">'Robust LS'</span>);

h = bar(hist3,freq3);
set(h,<span class="string">'FaceColor'</span>,<span class="string">'none'</span>);
h = bar(hist1,freq1);
set(h,<span class="string">'FaceColor'</span>,<span class="string">'none'</span>);

xlabel(<span class="string">'||(A0 + u1*A1 + u2*A2)*x - b||_2'</span>)
ylabel(<span class="string">'Frequency'</span>)
hold <span class="string">off</span>
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="wcrobls__01.png" alt=""> 
</div>
</div>
</body>
</html>